arXiv: astro-ph/0112276v 1 12 Dec 2001 


Mon. Not. R. Astron. Soc. 000, [I| pp| (2001) Printed 1 February 2008 (MN LMjTjX style file vl.4) 


Harmonic analysis of cosmic microwave background data 
I: ring reductions and point-source catalogue 

F. van Leeuwen , 1,2 A. D. Challinor , 2 D. J. Mortlock , 1,2 M. A. J. Ashdown , 1 ’ 2 
M. P. Hobson , 2 A. N. Lasenby , 2 G. P. Efstathiou , 1 E. P. S. Shellard , 3 
D. Munshi , 1 V. Stolyarov 1,2 

1 Institute of Astronomy, Madingley Road, Cambridge CBS 0HA, UK 
2 Astrophysics Group, Cavendish Laboratory, Madingley Road, Cambridge CBS 0HE, UK 
3 DAMTP, Centre for Mathematical Sciences, Wilberforce Road, Cambridge CBS 0WA, UK 


Accepted for publication in MNRAS: 11 December 2001 


ABSTRACT 

We present a harmonic model for the data analysis of an all-sky cosmic microwave 
background survey, such as Planck , where the survey is obtained through ring-scans 
of the sky. In this model, resampling and pixelisation of the data are avoided. The 
spherical transforms of the sky at each frequency, in total intensity and polarization, 
as well as the bright-point-source catalogue, are derived directly from the data reduced 
onto the rings. Formal errors and the most significant correlation coefficients for the 
spherical transforms of the frequency maps are preserved. A clean and transparent 
path from the original samplings in the time domain to the final scientific products 
is thus obtained. The data analysis is largely based on Fourier analysis of rings; the 
positional stability of the instrument’s spin axis during these scans is a requirement 
for the data model and is investigated here for the Planck satellite. Brighter point 
sources are recognised and extracted as part of the ring reductions and, on the basis 
of accumulated data, used to build the bright-point-source catalogue. The analysis of 
the rings is performed in an iterative loop, involving a range of geometric and detector 
response calibrations. The geometric calibrations are used to reconstruct the paths of 
the detectors over the sky during a scan and the phase offsets between scans of dif¬ 
ferent detectors; the response calibrations eliminate short and long term variations in 
detector response. Point-source information may allow reconstruction of the beam pro¬ 
file. The reconstructed spherical transforms of the sky in each frequency channel form 
the input to the subsequent analysis stages. Although the methods in this paper were 
developed with the data processing for the Planck satellite in mind, there are many 
aspects which have wider implementation possibilities, including the construction of 
real-space pixelised maps. 

Key words: cosmology: cosmic microwave background - methods: analytical - space 
vehicles - techniques: miscellaneous. 


1 INTRODUCTION 


The traditional method used to analyse the continuum dis¬ 
tribution of radiation on the sky is through pixelisation: the 
calibrated measurements are projected onto pixels of a sky 
map, and further analysis of the data is done by means of the 
pixelised responses. These methods, whicharemost power¬ 
fully explored in the HEALPix software ( |G6rski 1999 ), have 
the advantage of providing a visual link to the reductions. 
Another already well-developed example of pixelisation is 
the IGLOO scheme developed by Crittenden & Turok (1998) 
and Crittenden (2000). A properly designed pixelisation 


scheme can, in addition, simplify the further analysis, al¬ 
lowing the use of fast spherical transform algorithms, and 
removing contaminated regions of the sky (Gorski, 1994; 
Mortlock, Challinor & Hobson, 2001). A problem with any 
pixelisation of survey data is the loss of information due to 
the binning of the data. As a result of coverage variations, 
different pixels have different formal errors which compli¬ 
cates the application of fast spherical transforms. Pixels may 
be correlated with neighbouring pixels, which further com¬ 
plicates the analysis. The effective beam profile is unlikely 
to be circularly symmetric (as a result of intrinsic asymme¬ 
tries as well as due to unconvolved sampling intervals) and 
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requires the observations to be deconvolved. The effective 
two-dimensional beam profile applicable to a pixel therefore 
depends on the orientations of the scans which have con¬ 
tributed to it, information which is complicated to maintain. 
Thus, unless a range of ‘parallel maps’ with supplementary 
information is kept (the implementation of which would do 
away with most of the advantages of the pixelisation), a 
pixelisation of the data will inevitably lead to loss of infor¬ 
mation. This loss affects our ability to recover the statistical 
properties of the reduced data, and as such should be con¬ 
sidered a potentially serious disadvantage. 

The prod ucts of a survey mission like Planck (see 
f or example [h ttp://astro.estec.esa.nl/ planck ) or MAP 
( [http://map.gsfc.nasa.gov/m_mm.htm ) are a point source 
catalogue and a set of frequency maps, from which maps of 
astrophysical components are derived. Methods for the com¬ 
ponent separation have been developed amongst others by 
Hobson et al. (1998), Prunet et al. (2001) and Baccigalupi 
et al. (2000). In the latter paper an independent component 
analysis (ICA) algorithm is used to separate statistically in¬ 
dependent signals from the frequency maps. In Hobson et 
al. (1998) component separation on small rectangular fields 
is investigated using a maximum entropy technique in the 
Fourier domain. The latter method has been further devel¬ 
oped into a full sky imple mentation at the reso lution re¬ 
quired for the Planck survey ( |Stolyarov et al. 2002 ), with the 
frequency and component maps analysed in the form of their 
spherical transforms. The analysis of the cosmic microwave 
background (CMB; one of the components separated) de¬ 
pends ultimately on the reliability of the separation process. 
A major aim of the data analysis should be to derive the in¬ 
put for the component separation, the spherical multipoles 
of the frequency maps, in the most reliable way, preserving 
as much information as possible on formal errors and their 
correlations. 

The analysis of data from an all-sky survey mission in 
which the survey is built up from circular scans can be car¬ 
ried out efficiently in three steps, as was done for the Hippar- 
cos reductions (Lindegren, 1979; ESA, 1997): the first step 
reduces the raw measurements to ring data; the second step 
combines the ring data obtained over the mission to recon¬ 
struct the sphere; the third step analyses the sphere data. 
The Planck mission fits this model very well: data is obtained 
over 1 hour intervals during which the spin axis of the satel¬ 
lite remains in a nominally ‘fixed’ position. The collection of 
scientifically useful data starts at the end of a pointing ma¬ 
noeuvre and finishes with the start of the next manoeuvre, 
which moves the spin axis to the next scan position. During 
this interval, which is referred to as a time-ordered period 
(TOP), the satellite performs some 60 revolutions. The area 
of the sky covered by a detector during a TOP is referred 
to as a ‘ring’. The spinning of the satellite causes every de¬ 
tector d to describe a small circle on the sky, each with its 
own specific opening angle ad (the opening angle is equiva¬ 
lent to the co-latitude of the ring as seen from the spin-axis 
position). A slight wobble of the spin axis (nutation) widens 
the rings, but this effect is small relative to the beam width. 
Data in a ring are referred to the circle defined by the mean 
spin-axis position of the satellite during the TOP and the 
opening angle of the detector concerned. The nutation ef¬ 
fect will simulate a small-amplitude periodic variation of 
the opening angle for the detector, which will be most no¬ 


ticeable for point sources close to the rings as observed by 
the three highest frequency channels. Over a one-year mis¬ 
sion almost 9000 rings will be produced for each of the 48 
High Frequency Instrument (HFI) and 56 Low Frequency 
Instrument (LFI) detectors. 

In the harmonic model the three data analysis steps can 
be summarised as follows: 

(i) The analysis of the measured samplings per individual 
detector, collected over each TOP. The samples are referred 
to phases along reference circles of which the normal through 
the centre corresponds with the mean spin-axis position over 
the TOP. After cleaning these data from spikes (primar¬ 
ily resulting from very high energy radiation and particles), 
the phase-ordered data is collected in phase bins for fur¬ 
ther treatment. The short-term response variations are de¬ 
rived, using the accumulated differences between the mean 
responses in each phase bin and the individual contributions. 
After applying the response corrections, the mean sample 
values per bin are re-calculated and the data is ready for 
analysis. The first step in the analysis is the point-source 
transit identification. This is initially carried out on the 
phase-binned data, but in later iterations the point source 
catalogue constructed from the data is also used. The entire 
phase-binned signal is then fitted with Fourier components 
for the background signal and corrections to the assumed 
abscissae and intensities for all identified point sources. The 
corrected abscissae and intensities, with their formal errors, 
form the input to the point-source catalogue construction. 
The transit profiles of the brightest point sources are ac¬ 
cumulated to provide the input for the central-beam-profile 
calibration. 

(ii) Joint analysis of the Fourier components of the rings 
from all detectors operating at the same frequency, to pro¬ 
duce the spherical transforms of the frequency maps. For 
each ring and each detector there is a unique coupling ma¬ 
trix, which describes how the spherical multipoles are pro¬ 
jected on a small circle with a given opening angle ad and 
a two-dimensional beam profile. The accumulated coupling 
matrices and ring harmonics for a given frequency are the in¬ 
put to a generalised least squares solution for the multipoles 
of that frequency map. In an iterative process, the long-term 
response variations are calibrated; calibration of the outer 
beam profile may also be done at this stage. The processing 
up to this point is very similar for scalar and polarised data 
( |Challinor et al. 2002 ). 

(iii) The analysis of the spherical multipoles of the fre¬ 
quency maps (for example, by means of the maximum 
entropy method) to extract the faint point sources, the 
Sunyaev-Zel’dovich clusters (Sunyaev & Zeldovich, 1980) 
and the spherical multipoles of the component maps, which 
are further analysed in the form of their power spectra. 

Due to the way these steps are interlinked with geometric 
and response calibrations and the construction of the point 
source catalogue, they will require iterations before reduc¬ 
tion results can be considered satisfactory. 

The present paper focuses on the the first of the three 
steps mentioned above. Although this is seen here as the 
first step in the harmonic model for the data analysis, many 
considerations made in the reductions as described here are 
equally applicable in a more traditional pixelised approach. 
The second step is described in the accompanying paper 
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(Challinor et al. 2002), and covers the reductions of the 
ring data to spherical multipoles of frequency maps in de¬ 
tail, a process which is more specific to the harmonic data 
model. The third step, described by Stolyarov et al. (2002), 
is the full-sky maximum entropy component separation. Ad¬ 
ditional aspects concerning the analysis of incomplete sky 
maps are discussed in Mortlock et al. (2002). 

An important criterion for the implementation of the 
harmonic model is the stability of the spin-axis pointing. 
This stability is affected by two main contributors: residual 
velocities around the two axes perpendicular to the spin axis, 
which will create nutation; and external torques, which can 
create residual velocities. The amplitudes of the residual ve¬ 
locities at the start of a TOP are determined by the criteria 
for the nutation damping, and should be small enough not 
to create a significant effect (in comparison with the beam 
width) on the accumulated data. The nutation wobble of the 
spin axis and the effect of the external torques can be de¬ 
rived from an analysis of the satellite’s attitude in analytic 
form and through numerical integration. An analysis of this 
kind for the Planck satellite is presented in Appendix [X|, the 
results of which are summarised in Section LA 

In Section we show how the data analysis in the har¬ 
monic model is closely linked with geometric and response 
calibrations. Details of the data analysis method are pre¬ 
sented in Section Section |H] describes the construction of 
the point-source catalogue and the calibration of the refer¬ 
ence phases for the reduced TOPs. Section 5.2 presents the 
geometric calibrations: the focal plane geometry, the opening 
angle correction and the focal plane orientation correction. 
Finally Section describes the way the data reduction and 
calibration processes are linked in iterative loops. 


2 ATTITUDE ANALYSIS FOR THE PLANCK 

MISSION 

The results of a full analytic and numerical analysis of the 
Planck satellite attitude are presented in Appendix [\]. Here 
we summarize those aspects which are of immediate impor¬ 
tance for the data analysis. Though this analysis is applied 
to parameters associated with the Planck satellite, it can 
easily be modified to apply to different satellite configura¬ 
tions, as long as the outer product of the rotational velocity 
and angular momentum vectors is high with respect to the 
external torques acting on the satellite. 

The Planck satellite is designed for a survey mission: 
over a period of a year the entire sky will be scanned twice, 
providing maps of the microwave sky at pass bands with 
central frequencies between 30 GHz and 857 GHz. The scan¬ 
ning will take place at pre-determined nominal pointings of 
the spin axis (pointing in a roughly anti-Solar direction). 
These nominal positions will not be exactly reproduced: a 
pointing noise of a few arcmin is expected. The satellite will 
rotate at a nominal velocity of 1 rpm, but small variations 
are expected here too. The interval of ~ 1 hour between 
two re-positioning manoeuvres of the spin axis is referred 
to as a TOP. During a TOP, the satellite’s spin axis will 
describe a relatively stable nutation ellipse in the satellite 
coordinate system (for as long as it remains unaffected by 
nutation damping), the maximum amplitude of which is de¬ 
termined by the residual velocities around the x- and y-axes, 
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Figure 1 . The orientation angles in the Planck scan. The spin 
axis position (SA) is defined by its ecliptic coordinates (A 2 ,/3 2 ), 
measured from the vernal equinox (VE); the spin axis position 
and the North Ecliptic Pole (NEP) define a meridian, called the 
‘reference circle’; the position of a ring is defined by its opening 
angle ad] the scan phase ip defines the angle between the fiducial 
reference point (FRP) in the field of view (FOV) and the reference 
circle. 


and should be less than 1.5 arcmin at the end of the nuta¬ 
tion damping. The period of the nutation will be around 3 to 
4 minutes, and will slowly change over the mission due to 
changes in the inertia tensor. The noise on the movement of 
the spin axis of the satellite relative to the nutation ellipse 
is at the arcsecond level. Such variations are not relevant 
for detectors with a 5 arcmin or more fwhm response beam. 
This allows one to incorporate a simple set of nutation el¬ 
lipse parameters in the data analysis if necessary. 

The rotation of the satellite will cause every detector 
d to describe a small circle on the sky, at an angular sep¬ 
aration aa (referred to as the opening angle) from the po¬ 
sition of the spin axis (see Fig. nl). As the actual spin axis 
will not coincide exactly with the nominal spin axis as de¬ 
fined in the hardware, there will be a difference between the 
actual opening angle and its hardware specification value. 
The wobble of the spin axis slightly widens the ring, but 
for most detectors the beam width is such that this can be 
ignored. This may not be the case for the highest frequency 
detectors, where the spin axis wobble may add a little to the 
beam width as observed perpendicular to the scan direction. 
The way it does so will be unique to each ring set due to 
the variable conditions at the end of nutation damping (see 
Appendix [a]) . 

We can conclude that, within the model presented, the 
attitude of the satellite can be described in a simple manner, 
while the actual attitude motions do not cause any problems 
for the harmonic analysis of the data. The chance of signifi¬ 
cant unaccounted forces acting on the satellite, which could 
disturb this situation, is small. 
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3 THE HARMONIC DATA MODEL AND 
SYSTEM CALIBRATIONS 

The harmonic data model aims to provide an accurate anal¬ 
ysis of the data and its statistical properties. It does so 
through a sequence of reduction steps, resulting in frequency 
maps in the form of spherical multipoles with their covari¬ 
ance matrices and a list of point source positions and inten¬ 
sities with associated formal errors. These maps can then be 
used as input for the component separation pr ocess, which 
can be p erformed in spherical multipole space (Stolyarov et 


al. 2002 ). It should be possible to derive the spherical trans¬ 
forms of the frequency maps from the Fourier components 
of the rings, using information on pointing, focal plane ge¬ 
ometry, focal plane orientation, the opening angle, the noise 
correlations of the data and the beam profile, although a 
number of computational pr oblems makes this stage non¬ 
trivial (Challinor et al. 2002). 

The Planck mission is designed to obtain full-sky sur¬ 
veys at nine different wavelengths. It explores new ranges of 
resolution, coverage and sensitivity at these wavelengths. For 
such pioneering observations, calibrations are best done by 
demanding internal consistency of the observed sky, which 
in most cases will not have changed from one observation 
to the next. The exceptional variable and moving sources 
should all be identifiable. Demanding internal consistency 
defines simple, differential reduction procedures, and should 
result in mission products with the highest possible inter¬ 
nal precision. Absolute calibrations are carried out only in 
the final stages of the data processing, using a selection 
of sources that can be assigned and measured in a well- 
controlled way with ground-based or other suitable exper¬ 
iments. It is important during the absolute calibrations to 
consider the spectral response of each source and its convo¬ 
lution with the Planck passbands, and any background con¬ 
tributions which may be differently perceived by different 
instruments. The latter applies in particular to calibration 
sources close to the galactic plane. 

The preference for an internal, differential calibration 
applies to both responses and geometric properties. For the 
responses it involves the following steps: 

(i) The short-time response variations (time scales from 
one minute to one hour) are calibrated as part of the rin g 
reductions, using the phase-binned data (see Section 4.4). 
They show up by displaying, as a function of time, the differ¬ 
ences between the mean samples per bin and the individual 
contributions to each bin. Removing any systematic effects 
observed in these residuals reduces all data contained in one 
TOP to an arbitrary mean response level for a specific de¬ 
tector and ring. All variations at shorter time-intervals are 
reduced to noise contributions, some of which may be cor¬ 
related between neighbouring samples. This is taken into 
account in the noise correlation matrix for the model-fit so¬ 
lutions. 

(ii) The long-term response variations (with time scales 
from one hour to one year) are calibrated as part of the 
sphere reductions. The reconstructed spherical transforms 
will, in the first instance, represent an average of all mean 
ring responses. By relating (in Fourier space) the response of 
each ring to this mean sky, response corrections per ring and 
detector are obtain ed. This is the equivalent of de-striping 
( Rcvcnu et al. 200C ), though in the method proposed here all 


data on the ring will contribute. In case the instrument high- 
pass filters the signal, and there is essentially no information 
left on the monopole, these adjustments can still be done on 
the remaining signal. The response-corrected rings can be 
used again in a reconstruction of the spherical transforms, 
which should result in a decrease in the fitting noise. The 
process can be repeated until no further significant changes 
to the ring responses are observed. Given the expected large 
number of rings involved (9000 per detector), this process 
is expected to converge quite rapidly. The Hipparcos reduc¬ 
tions have demonstrated the validity of this approach in the 
iteration between the great circle reduction, the sphere re¬ 
duction and the improvements to the astrometric parame¬ 
ters, which was based entirely on internal consistency of the 
data (see for example van Leeuwen, 1997). External infor¬ 
mation may be used to obtain a reasonable first estimate of 
the calibration parameters to ensure a rapid convergence of 
the process. 

(iii) The last step in the response calibration adjusts each 
final spherical transform of a frequency map using a set of 
ground-based (or other) calibration sources (see above), thus 
transforming them from relative to absolute scales. 


The result of these processes is a set of frequency maps 
whose precision is entirely determined by the internal cali¬ 
brations, and accuracy by the quality and quantity of exter¬ 
nal calibration sources. At any stage, should improvements 
in external sources become available, the external calibra¬ 
tion can be performed again. In addition, as no external 
information is used in the internal calibrations, the loss of 
information from the original data is minimal. The loss is 
determined primarily by the limitations in the modelling of 
detector response variations. If instead, for example, prop¬ 
erties of the dipole are used in the earlier calibration steps, 
the extraction of independent information on the dipole from 
the Planck data could be compromised. The accuracy of the 
calibration parameters would depend on the local ampli¬ 
tude of the dipole signal as projected on a ring, which varies 
strongly with the ecliptic longitude of the spin axis. 

The response or beam profile forms the link between 
the detector response and the detector pointing. The cal¬ 
ibration of the central part of the beam profile is closely 
linked with the geometric calibrations and the construction 
of the bright point source catalogue. The geometric calibra¬ 
tions determine the pointing of each detector at any time 
during the observations. The relation between pointing and 
time can depend on the following elements: 


(i) The ecliptic coordinates of the mean position of the 
spin axis during a TOP: (A Z ,/3 Z ) (see Fig. |lj). This position 
and the position of the North ecliptic pole define a great 
circle, referred to as the ‘reference circle’. 

(ii) For each TOP the nutation ellipse parameters to, u>x o, 
uj z , fi and fi as defined in equation (B3) in Appendix |b| 
The first two of these determine the ^-amplitude and phase 
of the nutation ellipse, the other three determine the period 
and the ratio between the x and y amplitudes. The nutation 
ellipse is described in the satellite coordinates of the inertial 
reference system, as defined in Appendix When observed 
on the sky, the ellipse is itself rotating at the spin velocity of 
the satellite. The nutation period changes slowly due to vari¬ 
ation in /i and resulting from depletion of consumables 
on-board the spacecraft. 


© 2001 RAS, MNRAS 000, id 








Harmonic analysis of CMB data 5 


(iii) The spin velocity cj z and its linear dependence on 
time during a TOP. The latter is again a parameter which 
changes very slowly over the mission. 

(iv) For each TOP, the time of the first transit of the 
fiducial reference point (FRP) in the focal plane through 
the reference circle. A first approximation can be obtained 
from, for example, the star mapper data and the calibra¬ 
tion of the star mapper^] position with respect to the focal 
plane geometry. Corrections to the time of first transit are 
obtained during construction of the point-source catalogue. 

(v) The angle between the FRP and the mean spin axis 
position (the opening angle «o, where the index ‘0’ indicates 
the FRP rather than a specific detector ‘d’). This is a correc¬ 
tion to the ground-based calibration value, that will change 
slowly with time due to inertia tensor changes. 

(vi) The focal plane geometry, which describes, for an as¬ 
sumed orientation of the focal plane assembly, the coordi¬ 
nates of all detectors as projected on the sky, relative to the 
projection of the FRP in the focal plane. 

(vii) The focal plane orientation correction, describing 
the difference between the actual and the assumed orien¬ 
tation of the focal plane as a function of time. 

While items (i), (ii) and (iii) are probably best derived from 
the star mapper data, items (iv) to (vii) rely entirely on 
information contained in the science data, in particular on 
transits of bright point sources. Starting values for the focal 
plane geometry will be obtained from the instrument de¬ 
sign specifications, but still need to be verified in flight. The 
calibration of these parameters relies on observed abscissae 
and intensities of point sources for different detectors. The 
geometric calibrations are crucial to the reliability of the 
final mission products, and inevitably require a number of 
iterations through the data reduction and calibration pro¬ 
cesses. As a part of these iterations the bright point-source 
catalogue (BPSC from here on) is producec^. The BPSC 
provides the best possible relative positions and intensities 
of all brighter point sources detected by Planck, and is fi¬ 
nally used to calibrate the overall positional reference sys¬ 
tem to the International Celestial Reference System (ICRS; 
Kovalevsky et al. 1997). 

With the development of the BPSC, the calibration of 
the central beam profile can gradually be improved. In the 
first iteration through the data it is necessary to assume 
that all point sources pass through the centre of the beam. 
If, to a first approximation, the beam is represented as a 
two-dimensional Gaussian, then this assumption will do no 
harm: the width of the beam profile is in this case not a 
function of the ordinate of the transit. During subsequent 
iterations the ordinates of the transits can be calculated us¬ 
ing the geometric calibrations and point-source coordinates 
contained in the catalogue. Ordinate information can also 
be incorporated in the beam profile calibrations. The con¬ 
struction of the catalogue ensures that in the final iteration 

* A star mapper is a device recording the transits of bright stars 
for the purpose of spacecraft attitude control. 

T This is not the same as the early release point source cata¬ 
logue to be produced by Planck, which is intended to be an early 
compilation of approximate coordinates and intensities of point 
sources detected by Planck to be used in follow-up observations 
by different instruments 
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a complete sample of point sources (up to a maximum ordi¬ 
nate depending on the beam profile) are taken into account 
in the ring analysis. Th is is essential for the u se of such data 
in the sphere analysis JChallinor et al. 2002 ). 

The iterations through the data reductions are essen¬ 
tial due to the fact that many of the instrument calibration 
parameters have to be derived from the science data itself. 
This applies to both the response and the geometric cali¬ 
brations. The precision with which these parameters can be 
determined is an essential part in the final data-quality veri¬ 
fication. Poorly determined parameters will inevitably leave 
their mark on the final data products. The data analysis 
needs to identify any such parameters and their possible ef¬ 
fect on the data as part of its preparation for the scientific 
exploration. 


4 THE RING ANALYSIS 

The inputs to the ring analysis are the samples collected 
during a TOP for a single detector, supplemented by timing 
information (an absolute time for the first sample and a sam¬ 
pling length). It is assumed here that these intervals can be 
recognized from the satellite’s housekeeping data. The ref¬ 
erence position of the spin axis during a TOP is defined as 
the centre of the nutation ellipse actually described by the 
spin axis, and is derived from the star mapper data process¬ 
ing. Each detector d describes a small reference circle with 
opening angle ad around the reference spin axis position. 

The use of rings as an intermediate step in the data 
analysis of a full-sky survey mission was first proposed by 
Lindegren (1979) for the Hipparcos mission. The idea was 
independently explored for CMB surveys by Delabrouille, 
Gorski & Hivon (1998), who investigated the recovery of 
the CMB power spectrum directly from the ring data.. 

4.1 Signal representation 

The signal for the sky can be considered as consisting 
of three components: a continuous background, extended 
sources (which are not separated from the background) and 
point sources of various intensities. Translating this signal 
into the data observed on a given ring, a number of effects 
have to be taken into account: 

• The satellite is rotating at a scan velocity which is not a 
fixed value and which may in addition change slightly during 
a TOP; 

• The sky signal is convolved with a beam profile and 
sampled; 

• The sky signal has a (not necessarily constant) back¬ 
ground signal added to it, originating from the instrument 
itself; 

• The sky signal is represented on an arbitrary and pos¬ 
sibly slightly variable readout scale. 

Some of these effects are removed in the ring reductions, 
others in the co nstruction and calibrat ion of the harmonic 
frequency maps ( |Challinor et al. 2002 ). 

In the analysis of a TOP, two types of contributions are 
solved for: 

(i) The Fourier components, representing the structure 
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in the underlying continuum of the microwave sky and ex¬ 
tended sources: Co, C n and S„, with n = 1,..., n max ; 

(ii) The intensities and abscissae for point sources (or¬ 
dinates are assumed zero or obtained from the point-source 
catalogue): Ik and ipk, with k = 1,..., s, s being the number 
of sources solved for. 


In addition, the following properties of the data are resolved: 


(i) The detector response variations, such as changes in 
the background signal, A b(t — t T ), and drifts in the quantum 
efficiency, A q(t — t r ), where t r is an arbitrary reference time; 

(ii) The covariance matrix of the measurement noise. This 
will be the noise on b in-a veraged data rather than individual 
samples (see Section 4.4). 


The Fourier components for all TOPs corresponding to de¬ 
tectors in the same frequency channel are used to re construct 
the spherica l multipoles of the frequency maps (Challinor 
et al. 2002). This complements the “ring-torus” methods 


developed by Wandelt & Hansen (2001) for power spectrum 
estimation from the analysis of ring data from a special class 
of scans. The methods required for harmonic map-making 
also draw heavily on the harmonic-space convolution algo¬ 
rithms developed by Wandelt & Gorski (2001) and Challinor 
et al. (2000). The ring analysis is iterated with the construc¬ 
tion of the BPSC (see Section ^) bv means of the geometric 
calibrations mentioned in Section Q. The complete inclusion 
of identified point sources is also ensured through the BPSC. 
The success of this iteration determines the final pointing 
noise uncertainties and their contributions to the noise on 
the frequency maps. 


4.2 The time-to-scan-phase relation 


The first requirement for the processing of the TOP data is 
an accurate determination of the scan velocity and its change 
with time over the interval covered by the TOP. The scan 
velocity ui z {t ) determines the relation between time t and a 
relative scan phase %fi for the individual sampling intervals: 


tp(t) = ip(t r ) + 


u> z (t)dt. 


(1) 


The attitude simulations (Appendix |Xj) show that the scan 
velocity can to high accuracy be approximated by the scan 
velocity at reference time £ r , t o z {t T ) , and a constant scan- 
velocity drift, Cj z , as 


Wz(t) = Wz(tr) + {t~ t r )Wz, 


( 2 ) 


giving the following relation for ip(t)\ 

1p{t) = + (t — tr)0J z {tr) + i(t - tr) 2 W z . (3) 

A further disturbance on the time to scan phase relation 
comes from the wobble of the spin axis, but this is very 
small as is shown in Appendix |§J. 

The scan velocity and its variation will most proba bly 
be derived from the star mapper data. However, Section 4.8 
shows how information on the scan velocity is, in principle, 
also present in the science data. There may occasionally be a 
discontinuity in the relation between time and scan velocity 
due to semi-discrete torques caused by the satellite being hit 
by a micrometeorite. 


Using equation (|j), phases can be assigned to each sam¬ 
ple within the TOP. Sorting the data explicitly or implicitly 
(through a reference index) on ijj modulo 2-7T produces the 
phase-ordered data (POD) for a TOP. The use of POD has 
also been explored by Wandelt & Hansen (2001) in the con¬ 
text of power spectrum estimation. The further processing 
of the POD involves the following steps, which are described 
in detail in the sections below: 


(i) Spike detection is performed on the POD, where spikes 
will show up more clearly du e to the ~ 55 times higher 
density of data points (Section 4.3); 

(ii) Phase binning of the data: This achieves a very sig¬ 
nificant compression of t he data without significant loss of 
information (Section 4.4); 

(iii) Corrections for short-term detector response varia¬ 
tions (Section [l.5| ); 

(iv) Point source identification: either from the data 
stream itself or from the point-source catalogue, providing 
abscissae, ordinates and intensities (Section |l.8| ); 

(v) Signal fitting: solving for the Fourier components 
representing the continuum and corrections to some of 
the poin t-source parameters, and the noise spectrum (Sec¬ 
tion 4.9); 


Processes (i), (ii) and (iii) are applied only once to the data, 
while processes (iv) and (v) are part of the iteration with 
the BPSC construction. 


4.3 Spike detection and removal 

In the POD spikes will be much more conspicuous than in 
the TOD, as data points are compared with others at al¬ 
most identical telescope pointings. Filters have to be devel¬ 
oped that can reliably detect spikes as outlying points. All 
spikes are to be removed, while their times and intensities 
should be collected to allow for tests of their statistical prop¬ 
erties (distribution over time and intensity). It may also be 
necessary to remove the sample(s) immediately following a 
spike in the TOD, and to compare observed response vari¬ 
ations for detectors with the occurrence of spikes. Confu¬ 
sion between spikes and transient objects should be carefully 
avoided. Once the POD has been searched for, and cleaned 
from, spikes, it is ready for further analysis. 


4.4 Phase binning of the ring data 


A typical TOP for an HFI detector will contain around 
7 x 10 s samplings (200 Hz sampling frequency). The Fourier 
analysis is likely to require an £ max value of around 2500, 


giving some 5000 unknowns (see Section 4.7). The process 
of phase binning compresses the ring data without signifi¬ 
cantly modifying it (as would be the case when resampling), 
and as a result does not inflict a significant loss of infor¬ 
mation. The compression of the data, typically by a factor 
of 40 to 50, significantly reduces the processing time and 
data storage requirements. Relative to the original observa¬ 
tions, the phase-binned data provides a much improved ba¬ 
sis for the detection of point sources and detector response 
variations. The phase binning works for both point sources 
and the harmonic analysis of the background signal and is 
based on principles developed for, and used extensively in, 
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the Hipparcos data analysis of the 1200 Hz modulated main 
detector signal (ESA, 1997; van Leeuwen, 1997). 

The basic relation between an observation (one sam¬ 
pling Oi by a single detector) and its Fourier representation 
is given by 

n max 

Oi = Co + [ c " cos nipi + S n sin nipi ] + Ni , (4) 

71 = 1 

where ipi is the phase of the observation and Ni represents 
the instrument noise. We will use this representation to de¬ 
scribe the signal with all bright point sources removed. The 
ring is divided into m phase bins of equal length 27r /m. The 
phase at the centre of each bin is given by T,- = 2nj/m. 
Every observation is associated with a bin j, which turns 
equation (]!]) into 

n max 

Oi = Co + [C" cos n(4'_ ) + dipij) 

71 = 1 

+S„smn('t>j + dipij)] + Ni, (5) 

where dipij = ipi — T . This can be further expanded to 


Oi = C 0 +J2 


C n (cosn'l l j cos ndipi■ 


— sin n'bj sin ndipij) + S „(sin n't!j cos ndipi 


+ cos n4/ j sin ndipij ) 


+ Ni. 


( 6 ) 


Phase binning provides weighted means per bin j for the 
left- and right-hand sides of equation ([]). The weights for 
all contributions in a bin can be taken equal. If nj denotes 
the number of samples in bin j , then the following relation 
is obtained for the mean signal in a bin: 

. 71 max ‘ 

= -E°‘ = C7 ° + E 

i£j 7i=l L 

C„( E c , j cos n'hj — E s> j sinn^jJ 
+S n {f^a,j cos n't! j + E Ci j sinnTj) 

(7) 


+ JVj, 


where E Si) and E c j are defined below. By choosing the bin 
size sufficiently small, the sums over the bins can be es - 
timated. Experience with the Hipparcos data ( ESA 1981 ) 
showed that a coverage of a full cycle of the highest spatial 
frequency by 6 bins is sufficient, and allows for the following 
approximations to be made: 


Eg. 7 - 


Ec.7 - 


— ' sin ndipij 


iej 


1 


cos ndipij 


Y dipij = n(d'f/ j ), 
J iej 




1 - ^(d lp i: 


2 ny 


\2 i 2 2 

) = 1 - n a*. 


iej 


( 8 ) 


Depending on the amplitude of the nutation ellipse relative 
to the beam width, it may be necessary to fit the data con¬ 
tained in a phase bin as a second-order function of the offset 
from the reference ring position, primarily to avoid an in¬ 
tensity bias and noise increase due to point sources. 
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The phase binning thus requires the accumulation 
of four items per bin: the mean count Oj, the to¬ 
tal number of samples nj, the average phase correction 
<d*i> = £d</W nj and the phase dispersion ay ^ = 

y/y^( dipij) 2 /2nj. Values for (d4>) and ay are shown in 
Fig. Hfor a simulation using a sampling frequency of 200 Hz, 
n ma x = 2050, and 12500 bins. The offset from the nominal 
scan velocity was 1.243 arcsec s -1 , and the scan velocity 
drift was 0.009 arcsec s -2 . An average of 57 samples per bin 
is expected under these conditions. The expected value for 
(d^) is zero, and (ay) = ('K/m)/\/6 (where m is the total 
number of bins), which equals 21'.'6, equivalent to 5.8 per 
cent of the beam dispersion, or 5 per cent of the bin width 
in this example. 

The phase-binned data Oj and (dTj), ay^ and the cor¬ 
relation matrix of the A fj enter our model of the response 
through equations (Q) and Q in a generalised least squares 
solution for the {C n , Sn}- 

In the evaluation of equation (|s|) we ignored the fourth 
order term, which has an expected value of (n7r/m) 4 /120. 
With m = 6 n for the highest n value used, this amounts to a 
maximum correction of 6 x 10 ~ 4 to E Ci J . This figure should 
be compared with the expected amplitudes for the highest 
frequency harmonics (which will be severely depressed by 
the beam convolution), relative to the noise on the measure¬ 
ment Oj. For the lower frequencies the contribution of the 
fourth order term is very much smaller still (e.g. 6 x 10 -12 
for n = 25). In fact, for most of the lower frequencies the 
approximations of E s « 0 and E c as 1 can be used. 

Thus, phase binning brings down the number of obser¬ 
vations used for the ring analysis by a factor ~ 57, without 
any significant loss of information. The data storage require¬ 
ments are as a result significantly reduced. All trigonometric 
coefficients in the analysis are fixed in the binned solution, 
and can be pre-calculated. The phase-binned data for each 
ring can be kept as intermediate data products, as the con¬ 
tents do not change in the iterative processes. 


4.5 Short-term response variations 

As stated before (Section |4.l[ ), the short-term response vari¬ 
ations can be derived from the systematic differences be¬ 
tween the individual sample counts and the mean count of 
the phase bin to which the sample has been assigned. Cor¬ 
related noise between neighbouring samples will transfer to 
similar noise between neighbouring bins. Two types of cor¬ 
rections can be made: 

• A background variation correction A b(t — t T ), which can 
represent effects like changes in the radiation received from 
the instrument and the spacecraft, and which acts like a 
varying constant offset; 

• A quantum efficiency variation correction Ag(t — t r ), 
which will produce residuals scaled by the mean bin count. 

For a sample Oi in phase bin j with mean count Oj, this 
means 

A Oij = Oi — Oj = Ab(t — t T ) + Aq(t — t T )Oj + Ni. (9) 

The functional representation for these two corrections has 
to be decided upon early in the mission, on the basis of the 
features shown in the data. A fit with a low-order spline 
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Figure 2. Results for a binning experiment at the highest Planck resolution: 5 arcmin, using 12500 bins. On the left are shown the 
actual values per phase bin as observed, on the right the histogram of the distribution of observed values, with an arrow indicating the 
expected mean value. Top graph: number of samples per bin; middle graph: phase dispersion per bin; bottom graph: mean phase offset 
per bin. 


function will in most cases take care of any real variations. 
At this stage it is unnecessary, and potentially even damag¬ 
ing, to reduce the detector responses to an absolute scale, in 
particular if those calibrations would include a dependence 
on the phase angle ip, for which the reference phase has not 
been accurately determined yet. 

The observed variations should be compared with tem¬ 
perature records for the spacecraft and payload and the oc¬ 
currences of spikes, and correlations may be used in the re¬ 
moval of any observed variations. As a result of this calibra¬ 
tion, the data collected during each TOP will be free from 
short-term detector variations. 


4.6 Corrections for the satellite’s motion 

The Planck satellite will be positioned in a Lissajous orbit 
around the L2 Lagrangian point of the Sun-Earth system. 
It will therefore describe an almost circular orbit about the 
Sun with a one year period and a radius of ~ 1.01 AU. The 
period for the Lissajous orbit relative to L2 is around 179 
days, and will have an amplitude of around 10 5 km. Thus, 
the orbital velocity of the satellite is dominated by its mo¬ 
tion around the Sun and will be approximately 30.3 km s^ 1 . 


This causes fractional spectral shifts of AA/A « 10~ 4 , which 
is equivalent to 9 per cent of the dipole signal in the CMB ra¬ 
diation. The Planck mission aims at detecting much smaller 
anisotropies in the CMB, and these effects are therefore a 
significant distortion of the signal. The effect will be oppo¬ 
site in the two half-year surveys, and will be most noticeable 
near the ecliptic plane. 

The data can be corrected for this effect iteratively with 
the production of the frequency maps. The frequency maps 
can be prepared to relatively low values of Z max ~ 500 to 
produce all-sky spectral index maps. The velocity vector of 
the satellite together with estimated maps of the spectral 
gradient can then provide corrections to the observed inten¬ 
sities for each ring: 

AI ss A^-- cos#, (10) 

a A c 

where 6 is the angle between the velocity vector of the ob¬ 
server (which has magnitude v) and the observation direc¬ 
tion, c is the speed of light in vacuum, dl/dX is the local 
spectral gradient, and A I the local intensity correction. This 
effect then has to be integrated over the spectral response 
of the beam profile to correct the actual observed signal. 
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The positional effect, generally referred to as aberra¬ 
tion, is, to first order in v/c, given by 

sin Ad = ( v/c ) sin 9, (11) 

where A 9 is the difference between the propagation direc¬ 
tion of the radiation in a stationary reference frame and the 
actual moving reference frame. For the Planck observations 
A 9 has a maximum of 20 arcsec (for the part of the scan 
nearest to either of the ecliptic poles) and can in principle 
be taken into account when assigning scan phases to the in¬ 
dividual samples. This is likely to be relevant for the HFI 
detectors, for which the maximum correction is comparable 
with the abscissa accuracies for bright point sources, and the 
angular scale of the highest harmonics used in the signal 
analysis. Ignoring the correction will result in information 
leakage between neighbouring frequencies in the harmonic 
analysis, and a significant positional noise contribution for 
point sources near to the ecliptic poles. 

4.7 Resolution requirements 

An important aspect of the realization of the harmonic data 
model is the value of n max that should be applied to the 
data analysis. This value is largely determined by the beam 
width, the scan density, and the way in which point sources 
are dealt with in the analysis. Good estimates for n max 
can be obtained by assuming a circular-symmetric Gaus¬ 
sian beam profile. The power spectrum of a point source 
convolved with a Gaussian beam is shown in Fig. []. It is 
clear that in order to represent the point sources adequately 
as components in the harmonic analysis relatively high val¬ 
ues of n max are required. An alternative approach is to treat 
identifiable point sources as separate components in the ring 
analysis. This should be feasible if no features are expected 
in the background signal, which are sharp compared to the 
beam profile. 

A circular symmetric Gaussian beam with a dispersion 
(Tb (in radians) will reproduce the higher harmonics in the 
background signal with decreased amplitudes. If the beam 
is represented by 

= exp[--i/> 2 /2of], (12) 

then the decrease in the amplitude is given by 

aim' = aim exp[-(la b ) 2 /2], (13) 

where ai is the actual and aim the observed amplitude (see 
also Challinor et al., 2002). Provided that point sources are 
treated as separate objects, the suppression of higher har¬ 
monics determines the value for n max in the ring solution. 
For ad « 85°, simulations (to be detailed in a future pa¬ 
per) suggest a value of n max « 4400 aremin/at, where at is 
expressed in aremin. Thus, for a b = 2.12 aremin (fwhm of 
5 aremin) we find n max ~ 2100. At such high n max values 
the matrices involved in the transformations are very large, 
and the covariance matrix of the spherical multipole solution 
will contain around 10 13 elements. Fortunately, simulations, 
which will be detailed in a future paper, have shown that for 
many plausible scanning strategies, the covariance matrix of 
this solution will be very sparse (some analytic approxima¬ 
tions to the covariance matrices for simple scan strategies 
are also derived in Challinor et ah, 2001). 

© 2001 RAS, MNRAS 000, 





Figure 3. Bottom graph: The power spectrum of a point source 
(with sampled peak intensity equal to 1) given a Gaussian beam 
with FWHM of 5 aremin (a b = 2.12 aremin) and a sample width of 
1.8 aremin. Middle graph: the restored peak intensity of a point 
source if fitted with harmonic components up to the indicated 
n max value. Top graph: the restored point source for n ma x = 2500 
(thick), 4000 (intermediate) and 6000 (thin) (also distinguished 
by increasing peak height and decreasing side-lobe amplitudes). 


4.8 The point-source parameters 

The main tasks of the point source processing are the fol¬ 
lowing: 

• To identify from either the data stream or the BPSC 
the point sources present in the ring data; the source of 
information will depend on the iteration stage of the data 
reduction process; 

• To supplement this list with solar system objects that 
may have been observed; 

• To produce for each point source preliminary estimates 
of the intensity I s , the abscissa ?/> T , and if obtained from a 
catalogue or ephemerides (solar system objects), the ordi¬ 
nate tv; 

• To obtain as part of the reduction of the ring data the 
intensities I„ of all identified sources, and abscissae %p T for 
sources with a sufficiently high signal-to-noise ratio; 

• To collect the abscissa data for (re-)building the BPSC; 

• To collect the profiles and fitting parameters for the 
brightest sources as a contribution to the beam profile cali¬ 
bration. 

This process is clearly iterative, improving at each step the 
quality of the BPSC, the beam profile and the reliability of 
the segregation of the point sources from the data. 
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Figure 4. The point-source position (PS) relative to the phase 
binning and the mean scan circle. 


4-8.1 Identification of point sources 


The mechanism of point source identification will depend 
on the iteration phase of the data reductions: in the first 
instance point sources will mainly be identified from the 
data stream itself, while in subsequent iterations the iden¬ 
tification will rely more on the BPSC. It is expected that, 
at least for the high-frequency detectors, the power will be 
dominated at high l values by point sources, which should 
allow for the development of a reliable filter for the detec¬ 
tion of brighter point sources. Filters for the recognition of 
point sources in a one-dimensional data stream were devel¬ 
oped and success fully applied for th e Hipparcos and Tycho 
data reductions (van Leeuwen 1997). Simplified versions of 
two-dimensional filters under development for recognition of 
point sources on maps (see for example Cayon et ah, 2000, 
Sanz et ah, 2001) could also be considered. The detection 
of point sources from the data is significantly enhanced by 
the phase binning of the data, increasing the signal-to-noise 
rati o by almost a factor 7 for the parameters used in Sec¬ 
tion 4.4. 


On the first pass through the data all point sources are 
assumed to transit through the centre of the beam. This 
is not problematic if the beam is approximately circularly 
symmetric and Gaussian. The deviation of the actual beam 
profile from this assumption will cause a slight error in the 
first reconstruction of the beam. This error can be reduced 
once ordinate information on point source transits becomes 
available too. 

All detected point sources will be used to build the first 
and subsequent versions of the BPSC: the evolving catalogue 
used in iterations to predict and consistently identify point 
sources. In these iterations information on the ordinate of 
the source at the time of the transit should also be incorpo¬ 
rated. The consistent inclusion of point sources in the ring 
analysis is an essential requirement for the further analysis 
of the underlying continuum. The BPSC also plays a crucial 
role in the geometric calibrations of the instrument. 


4-8.2 The measured and binned point-source signal 

The sampled responses Oi for a point source of intensity I s , 
passing through the beam of a detector, is a function of the 
scan phase -)/>; (at the midpoint of the sampling interval) for 
sample i and the abscissa ip T and ordinate v r of the point 
source (see Fig. 0): 



i p (arcmin) 

Figure 5. The beam response R(tp, 0) (thin line) and the con¬ 
volved beam response S(ii). 0) (thick line) profiles for a Gaussian 
beam with FWHM= 5 arcmin. The increase in beam width as a 
result of the samplingis of the order of 20 arcsec. Top graph: in 
linear response scale; bottom graph: in logarithmic response scale. 


/“i/jj+A-j/j 

Oi = I s / - ip,v T )dip, 

J 'ijji —A ijj 


(14) 


where the integral represents the sampling interval, and 
R(i/Jt — ip, Vt) is the normalized beam profile for a detec¬ 
tor as a function of the offset from the centre of the beam 
(ip T —ip along the scan direction and v T perpendicular to the 
scan direction). The effective beam profile S(ip,v) is defined 
as the actual beam profile R(ip, v) convolved with the sam¬ 
pling interval, as illustrated in Fig. 0 [see also equation (0)] 


(15) 


S{ip, v) = / R(ip',v)dip', 

J — 

where the relevant coefficients are defined in Section 4.4 
When phase binning is applied to the point source contri¬ 
butions, the response in phase bin j is given by 


Oj = I s [S(lpT — ^j,V T ) — (d'L 3 ')5’ , (^>r — 'S/jjVr) 

+al.S"(iP T - ® 3 -,u T )] +Bj + Mj, (16) 


where Bj represents the background signal. In the same way 
as was found for the harmonic multipoles, the main contri¬ 
bution comes from the second derivative, which produces an 
effective broadening of the beam. However, if we use 12 500 
bins for a beam with fwhm=5 arcmin (at, = 2.12 arcmin), 
the additional dispersion is u^ j « 21” 6 and the effective 
fwhm of the beam is increased as a result of the phase bin¬ 
ning by no more than one per cent. 


4-8.3 Fitting parameters for the point-source signal 

Two parameters require determination in the signal fit: the 
transit phase ip T and intensity I a of the point source. For 
both, preliminary estimates are required which are then ad¬ 
justed in the solution. 

Equation [1^ can be linearized in the intensity correction 
SI s and the transit correction Stp-r so that 

Oj = SI s S(ip T — ^ j,v T ) + Sip T IsS'(ip T — 'Lj, «r) 

— 1 -I s S(ip T — 1 Fj, v T ) + Bj + Afj. (17) 
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This equation enters in the ge nera lised least squares solution 
for the ring data (see Section 4.9). The intensity correction 
term in Eq. [ujis scaled by the response function. As a result 
of this, the noise on the lower intensities (the wings of the 
response function) will tend to affect the determination of 
5I S more than the better determined higher intensities (core 
of the response function). This is in particular the case when 
the noise on the signal depends on its intensity. Therefore, 
only a few of the central phase bins should be used for solv¬ 
ing for the transit parameters of a given source. It may be 
necessary to iterate the ring solution to properly separate 
the point source and continuum contributions. 

Assuming we use the central 5 phase bins for determin¬ 
ing the fitting parameters for each point source, estimates 
can be obtained for the expected precisions. For the inten¬ 
sity error a value of 0.8 times the noise level on the binned 
samples is expected, equivalent to 0.12 times the noise level 
on individual samples. For the abscissa error we estimate a 
value of 1.7 arcmin divided by the signal-to-noise ratio of 
the point source in the phase-binned data. 


4-8.4 Transient sources 

Any object moving by a significant fraction of the beam 
width over an interval of one hour will have to be treated 
as a transient source. For the highest frequency channels 
of Planck this translates into a movement of approximately 
2 arcmin per hour and above. The fastest objects in the 
Earth’s neighbourhood will be traveling at ~ 10 5 km hr^ 1 . 
The limit of 2 arcmin then translates into a horizon of 
around 1.2 AU, which will include a small fraction of as¬ 
teroids and the occasional comet. 

For sources that are not variable on a time-scale of less 
than one hour (which may not be true due to rotation of 
the objects), systematic differences between the actual mea¬ 
surements (before phase binning) and the reference profile 
can be expressed as corrections to the effective scan velocity 

0Jz{tr) ■ 

AOi = Oi — Bi — I s S{ipT —ipi,v) 

F) Q 

= {t-T)I s -^du z + Ni, (18) 

which shows that most of the information on u} z (t T ) is con¬ 
tained in those sections of the signal which have the steepest 
gradient as a function of ip. The same kind of information 
can in principle also be used to determine the scan veloc¬ 
ity from the science data using bright-point-source measure¬ 
ments, although it would be by far preferred to derive this 
information from the star mapper data. 


4-8.5 Beam profile calibration 

The beam profile calibration uses the accumulated transit 
data of bright point sources (which may include transits of 
solar system objects). During the first processing of the data 
no reliable ordinate information is available (catalogue posi¬ 
tions are still poor or non-existent, and opening angles still 
need to be calibrated), and every source is assumed to go 
through the centre of the beam. During later stages of the 

© 2001 RAS, MNRAS 000, 


data reductions, when the point-source catalogue and the 
geometric calibrations are improving, the ordinate informa¬ 
tion can be incorporated too. 

The principle of the beam profile calibration follows 
techniques used in the Hipparcos data reductions for the 
single slit response functions of the star mapper slits (ESA, 
1997; van Leeuwen, 1997). The response function is sam¬ 
pled at a frequency 4 to 8 times higher than the sampling 
frequency of the data. Using the abscissae and intensities 
of identified point sources, the observed residual samplings 
Oj — Bj (original signal Oj minus the background Bj as es¬ 
timated from the current harmonic fit to the continuum) in 
phase bins close to the point source transit phases {ip T } are 
binned according to their separation from the transit phase 
{ip T }- Also accumulated in bins are the derived transit peak 
intensities I s . Thus, we find for the beam-profile bin k the 
mean normalised response 


S(k) 


£ [©,■(*)-»,-(*)] 

E B(k) 


(19) 


where the sums are over the point sources and phase bins 
contributing to the beam-profile bin k. Noise correlations 
between neighbouring bins in the phase-sampled data will 
also affect the accumulation of equation (|lR]), in that pairs 
of bins in the accumulation can contain partially correlated 
noise, which needs to be taken into account when fitting a 
response curve. The values S(k) can be fitted with a cubic 
spline to provide a smooth beam profile with a continuous 
derivative, which then provides an estimate for the sampled 
beam profile S(ip, 0). When, after the construction of the 
first point-source catalogue, ordinates become available for 
the point-source transits, the beam profile can be resolved 
in two dimensions. The reconstructed beam profile obtained 
this way is the sampled, effective profile, and not the actual 
profile, which would apply to a stationary detector. 

Transits of planets such as Jupiter may also be useful 
for the beam profile calibration, though could be problem¬ 
atic: the detector response to very high intensities will not 
be linear, and the spectral gradient for the planets is likely to 
be quite different from that of the average microwave point 
source. The beam profiles will vary with frequency (see for 
example Challinor et al., 2001), making it still more compli¬ 
cated to incorporate the profiles measured from the planets. 


4.9 The ring solution 

After binning the data, calibrating and removing short-term 
detector responses, and identifying the point sources, the 
TOP is ready for reduction. This part of the reductions 
consists of a generalised least squares solution for the ring 
harmonics in the underlying continuum (equations (pj) (^)) 
and simultaneous solutions for the abscissa and intensity 
corrections for all identified point sources (equation ([Tt])). 
The observation vector 2 has as its components the mean re¬ 
sponse in each phase bin. The observations are related to the 
vector x containing the amplitudes of the circle harmonics, 
{CVi, Sn }, and the corrections to the point source parame¬ 
ters, 5I s and Sip T , and the noise vector v, whose components 
are the noise in each phase bin, A/}, via the linear equation 


z = Aa; + v. 


( 20 ) 
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The matrix A depends on the locations of the phase bins, t&j, 
the phase corrections and dispersions, (dJ/-,) and a<a j , the 
current estimate of the point source intensities and positions, 
I s and ip T (and ordinate information as this becomes avail¬ 
able), and the current estimate of the beam profile, v). 
Equation ( po| ) is the matrix form of equations (jrj) and (pd|). 
The minimum-variance estimate of x is 


X = (A t N' 1 A)" 1 A t N“ 1 z, 

(21) 

with errors 


((x — x)(x — x) T ) = (A r IM _1 A) -1 , 

(22) 


where N = ( vv T ) is the (phase-binned) noise covariance 
matrix. 

Assuming the instrument noise is a stationary, random 
process with correlation function C(t ) in the time domain, 
we can obtain the noise contribution to the mean response 
in a given phase bin by integrating over the time periods 
corresponding to those observations falling in that bin. The 
covariance matrix [N] = (AfjAfji) then takes the form of 
a convolution: 

[NW = {jQ 2 Hx)C{T s [x + (j - j')/m]} dx, (23) 

where N s is the number of times the ring is scanned in 
one TOP, T s is the average spin period in that TOP, m 
is the number of phase bins, and the integration limits 
x± = —(N s — 1) ± 1/m. The function A(*) is given by 


A(x) = 


E 

n=-(JV s -l) 


0 (l/m — \x — n|) 


x (1/m — \x — n\)(N s — |n|) 


(24) 


where 0(a:) is the Heaviside unit step function, and arises 
from the effects of phase binning and repeatedly scanning 
the ring. For white noise the A/} are uncorrelated, but in 
the presence of a significant low frequency component cor¬ 
relations will arise. Typically, instruments are designed with 
the goa l of restricting coloured noise to frequencies below the 


spin frequency, in which case the correlation length of the 
noise exceeds the spin period. (For Planck HFI, the nominal 
knee frequency at which the character of the noise changes 
from 1// to white is less than 0.06 rad s -1 , while the spin 
frequency is O.lOrad s^ 1 .) As the correlation length becomes 
large compared to the spin period, the noise covariance ma¬ 
trix approaches Wn' = X 2 c + Xufijj', corresponding to a 
fully correlated offset in every bin with r.m.s. \c, and un¬ 
correlated noise with r.m.s. Xu- In practice, the noise power 
spectrum will have to be estimated from the data directly 
rather than relying on simple parameterised forms like that 
give above. If the correlation length exceeds N s T a the off¬ 
sets on different rings will generally also be correlated, so an 
optimal analysis would require reducing several rings simul¬ 
taneously. 

Under the assumption of Gaussianity, the relevant com¬ 
ponents of the covariance matrix ((x — x)(x — x) T ) determine 
the errors on the ring harmonics (marginalised over the point 
source corrections). We have conducted several simulations 
to investigate the effects of variations in the spin velocity, 
and the presence of point sources, on the covariance matrix 
for the errors on the ring harmonics. To isolate the effects 



Figure 6. A histogram of the normalised, absolute values of the 
off-diagonal elements in the square root of the covariance matrix 
for a typical bin-sampled ring solution, using l max = 512. The 
matrix can effectively be regarded as diagonal. 


of spin velocity and point sources we have only considered 
white instrument noise. For our simulations we adopted a 
maximum value of n equal to 512. The covariance matrix 
for the errors on the ring harmonics is found to be very 
close to diagonal if only relatively faint point sources are 
present. In the presence of a very bright point source (like 
Jupiter), which will create a very low-weight gap in the ring 
data, there is a minor effect on the diagonal structure of 
the covariance matrix, the effects being largest in low fre¬ 
quency detectors. The normalised off-diagonal elements in 
the (Cholesky) square root of a typical covariance matrix 
are accumulated in a histogram in Fig. || The results show 
that the correlations between the different ring harmonics 
are at a level of 0.2 per cent or less. These levels decrease 
still further with increased resolution. The inclusion of a 
bright point source would affect the distribution, but corre¬ 
lations would still be below the 1 per cent level. If we further 
include stationary instrument noise, we can expect the er¬ 
rors on the ring harmonics to still be uncorrelated between 
Fourier modes for large N s . If coloured noise is confined to 
frequencies well below the spin frequency, cor relations be- 
tween rings will only affect the n = 0 modes (Challinor et 
al. 2002|) . 


5 THE BRIGHT-POINT-SOURCE 
CATALOGUE 

The bright-point-source catalogue (BPSC) is both a major 
product of the Planck mission, and an important calibration 
tool. The catalogue is constructed iteratively from the ab¬ 
scissae and intensities of the point sources obtained in the 
current ring reductions. Corrections to the BPSC will ul¬ 
timately provide positions, intensities (at different frequen¬ 
cies), and variability information for all detected sources. As 
was explained in Section [], the spin axis position (including 
nutation movements) and spin rate are most likely derived 
from the star mapper data, all other geometric calibrations 
depend at least to some extent on the science data. 


5.1 The positional sphere reconstruction 

The construction of the BPSC follows a simpli fied versio n 
of the Hipparcos sphere reconstruction process (ESA 1997). 
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This process requires a set of a priori positions for the point 
sources, to which corrections are determined. These posi¬ 
tions are also required in order to identify consistently point- 
source transits on different rings with objects on the sky. 
Initial positions at a precision of twice the beam size (10 ar- 
cmin for Planck) will be sufficient for this purpose. Using 
the ground-based specifications of the focal plane geometry 
and inertia tensor, first estimates are obtained for the open¬ 
ing angle and focal plane rotation corrections, and the fo¬ 
cal plane geometry. These estimates can be further refined 
through the use of planetary transits. From this informa¬ 
tion we obtain preliminary details of the alignment of the 
different detectors during a scan: phase offsets and effective 
opening angles. The transits recorded in scans for the dif¬ 
ferent detectors can now be transferred to a common sky, 
where each transit is represented by an error ellipse, strongly 
elongated in a direction perpendicular to the scan direction. 
A simple algorithm is then required to cross identify the 
different transits, and to produce from the accumulated in¬ 
formation the necessary initial positions. 

Once a catalogue of initial positions has been obtained, 
and transits have been identified with sources in the cata¬ 
logue, corrections to the assumed positions and to the refer¬ 
ence phases of the rings can be determined. Given the eclip¬ 
tic coordinates (A z , f3 z ) of the spin axis, and a position for a 
point source i at (Xi,Pi), the expected abscissa and ordinate 
of the point source with respect to the ring can be derived. 
The coordinates ( ipi , pi) in the system with the spin axis at 
the pole (where ipi is measured from the reference circle as 
defined in Fig. [I], and Q is the latitude of the source relative 
to the great circle associated with the spin-axis position) are 
obtained from 


" COS Ipi COS (pi " 

= R 2 (/3 2 -J)R 3 (-A z ) 

' cos A i cos pti " 

sin ipi cos (pi 

sin Ai cos Pi 

sin pi 


sin Pi 


where R pep) is the matrix representing a right-handed ro¬ 
tation around axis i through angle <p. Equation ( pB| ) can be 
used to derive the relation between corrections to the source 
position (A;, f3i ) and the resulting changes to the scan phase 
ip (obtained from the transit abscissae) and ordinate v (re¬ 
flected in the transit intensities). We define p = ^ — a + v, 
where v is the ordinate relative to the actual ring, and a 
the opening angle. The relations between the offset values 
for (Aipi, Avi) and (AA;, Af3i) are then found to be 

A tpi cos 2 C* = fi AAi cos j3i + f 2 Aft, 

AwiCosCi = — f 2 A\i cos(3i + fiA(3i, (26) 

where /i and f 2 are defined as 
fi = sin (3 Z cos /3i — cos /3 Z sin Pi cos(A; — A z ), 
f 2 = cos/3 z sin(Ai — A z ). (27) 

For a scanning strategy where the spin axis remains in or 
close to the ecliptic plane (and Ai — X z ~ ±7t/2), /i « 0 
and I/ 21 ~ 1 for most sources, the exception being sources 
situated close to either of the ecliptic poles. The result is 
that very little information on the ecliptic longitudes can be 
extracted from the measured abscissae, except for images at 
high ecliptic latitudes. Some information on the longitudes 
can, however, be extracted using the measured intensities, 
which depend on v through the beam profile. This method 
will inevitably fail when the source is variable. 
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The positional sphere solution consists of a least squares 
fit of the differences between the observed and predicted ab¬ 
scissae, Atpij (for point source i on ring j) to corrections 
(AAi, A (3i) for each point source. In addition, there is a cor¬ 
rection Aipj to the assumed reference phase for each ring. 
In the solution a boundary condition needs to be included 
which states that the average correction ^ Aipj = 0 , else a 
singularity may occur. When using different detectors, the 
phase offsets between the detectors have to remain effec¬ 
tively fixed: the only changes can come from rotation varia¬ 
tions of the field of view, and focal length variations of the 
telescope, and both effects are unlikely to be significant. 

To estimate the number of observations and parameters, 
we assume 10 000 point sources on the sky bright enough 
to be detected on single rings (Barreiro, private communica¬ 
tion), and 10 000 rings per detector over the Planck mission. 
With the beam’s fwhm= 5 arcrnin, the average number of 
point sources detected per ring is 6 passing within ±<7t,, and 
6 more passing within ±2<T6 from the mean position of the 
ring. The total number of transits recorded is then ~ 10 5 
per detector, giving ~ 10 observations per source per de¬ 
tector. The total number of parameters to be estimated is 
~ 3 x 10 4 , and this does not increase when information from 
more detectors is used (except when some point sources are 
not visible for all detectors used). Increasing the number 
of detectors increases the number of observations and will 
increase the rigidity of the solution. 

The relative positions of point sources obtained this way 
are rigid, but absolute positions are only determined up to 
an overall rotation, and so requir e linking to the Interna - 
tional Celestial Reference Frame (Kovalevsky et al. 1997). 
This can be accomplished through cross identification with 
radio and optical counterparts, thus determining a global ro¬ 
tation to be applied to the entire catalogue. This may not be 
important for the statistical properties of the CMB compo¬ 
nent in the background, but is relevant for relating sources 
and the dust map to other observations. 


5.2 Geometrical calibrations 

In Section ^ we mentioned four geometric calibrations which 
rely on the point-source data collected by the satellite: 

(i) The reference phase of the FRP, to be obtained for 
each circle; 

(ii) The opening angle for the FRP, to be obtained as a 
slowly changing function of time; 

(iii) The focal plane orientation, to be obtained as a 
slowly changing function of time. 

(iv) The focal plane geometry, to be obtained as a fixed 
set of parameters. 

The first three points concern the first order geometric 
orientation of the field of view: its shifts perpendicular to, 
and along, the scan direction, and its rotation. These result 
in systematic shifts in the abscissae and blurred intensity 
distributions for point sources, both as functions of the field 
of view position of the detector. The systematic shifts can be 
solved for as instrument parameters in the positional sphere 
reconstruction. This also applies to the calibration of the 
along-scan position of each detector in the focal plane ge¬ 
ometry. 
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The opening angle for each detector can only be de¬ 
rived from the brightness distributions of the point sources 
observed with it. A maximum likelihood solution, optimiz¬ 
ing the intensity distributions of the brightest, non-variable 
point sources should provide these calibration values. This 
can, however, only be applied to sources for which the po¬ 
sitions on the sky are well determined in both coordinates 
through the positional sphere reconstruction. For the Planck 
mission this condition limits its application to point sources 
near to the ecliptic poles. 

Additional information can be obtained from transits 
of planets and minor planets, for which the absolute coordi¬ 
nates at any time during the mission are known to a much 
higher accuracy than required for the Planck calibrations, 
and which have the advantage of being visible to most or 
all detectors. To use these transits for opening angle cali¬ 
brations requires, however, accurate knowledge of the beam 
profile perpendicular to the scan direction as well as accu¬ 
rate predictions of their expected intensities, both of which 
may be difficult to obtain. 


Although the methods presented in the current paper 
were developed as part of the harmonic data model, most 
would also be useful when using the more traditional pixeli- 
sation methods. This applies, for example, to the iterative 
cycle of the ring reductions, BPSC construction and the ge¬ 
ometric calibrations. Phase binning of the data can also be 
used with pixel-based methods, as it provides a means for 
short-term response calibrations, point source recognition 
and data compression. 

Further work is in progress in areas of point-source 
recognition from the ring data, and the cross-identification 
of point sources as detected on different rings. 

The part of the Planck data processing presented in this 
paper will be very demanding. An estimated half a million 
rings will be produced by the HFI per year of observations. 
The calibration requirements will make it necessary to pro¬ 
cess each ring at least three times to get all geometric and 
beam profile calibrations implemented. The processing of 
such large quantities of data requires careful planning. 


6 ITERATIONS WITH CALIBRATIONS AND 

THE RING ANALYSIS 

Iterations between the point-source catalogue and the ring 
reductions are necessary to obtain a properly calibrated geo¬ 
metric reference system for the observations. Without these 
calibrations in place interpretation of the data in the form 
of (partial) maps will be of limited value, especially towards 
the higher frequencies in the power spectrum. 

Iteration with the BPSC is also required to assign ordi¬ 
nates to point-source transits, which is essential in calibrat¬ 
ing the two-dimensional beam profile. Using the BPSC for 
the identification of point sources in the final ring reductions 
ensures that the background signals for all all rings contain 
compatible information. Inconsistent point-source subtrac¬ 
tion would lead to harmonic signals in the ring analysis that 
can not be combined in the harmonic map-making; _ 


An itera tion with the harmonic map-making (Challinor 
et al. 2002) for the half-year data is required to remove 


the spectral shift in the data which results from the veloc¬ 
ity vector of the satellite (see Section [T(|. While the CMB 
dipole affects only the CMB component in the frequency 
maps, the satellite’s velocity vector affects all component on 
every frequency map in both aberration and Doppler shift. 

Most of these iterations require complete reprocessing of 


the 


With- 


nna data and mu, as such, very time consuming, 
out the m, however, the scientific interpretation of the. data 
will be subject to considerable uncertainty. 


7 CONCLUSIONS 

The harmonic data model, of which the first part was pre¬ 
sented in the current paper, provides a high level of infor¬ 
mation preservation in the data reductions. It defines cali¬ 
bration requirements and methods and provides a clear path 
from observed quantities to the scientific products. The lat¬ 
ter is important, as the interpretation of the science products 
requires reliable knowledge of their statistical characteris¬ 
tics, which are well defined in the harmonic model. 
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APPENDIX A: PLANCK ATTITUDE 
DYNAMICS 

The dynamics of the Planck satellite are largely determined 
by the following characteristics: 

• The total mass and its variation over the mission; 

• The inertia tensor and its variation over the mission; 

• The position of the centre of gravity (CoG); 

• The actual scan velocity, which will be very close to a 
fixed nominal value; 

• The size, position and optical characteristics of the solar 
panel. 


A1 Input parameters 

Early models by Matra Marconi (Dynamics and Pointing, 
260/CDA/NT/83.95) have provided an initial set of values 
for the satellite characteristics, which we can use to derive a 
model for the satellite attitude. Although changes in these 
parameters can be expected for the final realisation of the 
Planck satellite, the general character of the results pre¬ 
sented here will not be much affected. 

We use the following values: 


• The total mass is 772 kg; 

• The inertia tensor (defined in the SRS coordinate sys¬ 
tem, see the next section) is given by 



'699 

4.0 

4.5' 

= 

4.0 

766 

4.2 


.4.5 

4.2 

970. 


where the actual values used are only of relevance in as far 
as that they reproduce the approximate characteristics of 
the tensor; 

• The position of the CoG, in the same reference system 
as defined above, is given by 


"* 0 " 


' 0.0319 ' 

yo 

— 

-0.0315 

_Zo_ 


. 0.7282 . 


• The diameter of the solar panel is 4.5 m, with a secu¬ 
lar reflection coefficient C s = 0 . 17 , and a diffuse reflection 
coefficient Ca = 0 . 10 ; 

• The solar radiation pressure is 4.5 x 10 - 6 Nm -2 . 


The solar panel shields the rest of the satellite from 
solar radiation. Its position, size and optical characteristics 
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therefore determine the solar radiation force experienced by 
the satellite. The position of the CoG then determines how 
much of this force is experienced as a torque (see for example 
Spence, 1978). 


A2 The inertia tensor and coordinate alignments 


The inertia tensor given above is not diagonal, i.e. the prin¬ 
cipal axes of the satellite do not align with the principal axes 
of the inertia tensor. We define two reference systems; the 
Satellite Reference System (SRS), aligned with the princi¬ 
ple axes of the satellite, and the Inertia Reference System 
(IRS), aligned with the principle 2 -axis of the inertia tensor. 
The origin of both systems is the CoG. 

The SRS is defined as follows. The 2 -axis is normal to 
the solar panel, goes through the CoG and is pointing away 
from the Sun. The axaxis lies in the plane through the z- 
axis and the vector from the CoG to the FRP as projected 
onto the sky, and is orthogonal to the 2 -axis. The y -axis 
completes the right-handed triad. 

The IRS is defined such that the off-diagonal elements 
I xz and I yz of the inertia tensor are zero. This is obtained 
by two small rotations ai, 2 . A third rotation ( 03 ) is added 
to make z ', x' and the direction of the FRP on the sky 
coplanar. The three rotations are clockwise around the x, y 
and 2 axes respectively, and define a rotation matrix 


R « 


1 

Ct3 
—a 2 


— Oi3 

1 

Ql 


OL2 
— OLl 
1 


(Al) 


The inertia tensor I' in the IRS is related to the inertia tensor 
I in the SRS through 

I' = RIR t . (A2) 


The physical meaning of the angles 0:1,2,3 is simple: oi is 
the rotation of the focal plane relative to the ring; 02 is the 
difference between the nominal and actual opening angle for 
the centre of the focal plane; 03 is an overall phase shift. 

Due to the depletion of consumables, the values in the 
inertia tensor cannot be assumed constant. The angles 01,2,3 
will change over the mission as a result of this, 01 and 02 will 
require calibration during the mission, while 03 corrections 
are absorbed in the zero-phase corrections that have to be 
applied to each circle. 


A3 The attitude dynamics 

The attitude dynamics of a satellite describe the motion of 
the satellite around its CoG in a suitably defined reference 
system. We consider the satellite to be a rigid body, but the 
effects of deviations from this assumption still need to be 
investigated: a possible source of violation of the rigidity as¬ 
sumption is the circulation of cooler liquids. As a rigid body, 
the satellite’s motions are described by the Euler equation, 
defined here in the IRS: 

= JV'-W x I'u, (A3) 

where N' represents the external and internal torques acting 
on the satellite, and w the inertial rates around the axes of 
the IRS. Of the three components of u, is by far the 
dominant one, with a nominal value of 0.1047 rad s -1 . As 
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a result, the cross product in equation (A3) is only relevant 
for the x and y coordinates. This allows us to approximate 
equation (A3) as 


date 

df 

duty 

d t 
duj z 
dt 


K/I'xx - flUyUJ z + jr-UxUz, 

-‘■XX 

/' 

N ' y /I'yy + hu x UJ z - -j^-UJyUz, 

lyy 


N'z/I'zz 


(A4) 


where the third terms in the first two of these equations are 
ignored in the analytic model, and where 


fi 


h 



(AS) 


The quantity A = 1 + VJ 1 J 2 is generally used in specify¬ 
ing the dynamic characteristics of the inertia tensor. In the 
preliminary specifications for the Planck satellite A « 1.35. 
Here we use instead 7 = A — 1 in the development of the 
dynamic equations. 


A4 The solar radiation torque 


The solar radiation acts only on the circular solar panel, 
which is shielding the satellite. Relative to the Sun, the satel¬ 
lite’s orientation is described by the solar aspect angle £ be¬ 
tween the spin axis and the direction of the Sun, and the 
spin phase ip, measured in the direction of rotation from the 
transit of the focal plane through the great circle defined by 
the directions of the Sun and the spin axis (see Fig. 0). This 
defines the unit vectors s' in the IRS and s = R _1 s in the 
SRS, from the solar panel in the direction of the Sun as 


" — cos ip sin £" 


" — cos ip sin £ + 02" 

sin ip sin £ 

, s w 

sin ip sin £ — ai 

— cos £ 


— cos £ 


where £ is assumed to be very small (< 0.03 rad.). Following 
Spence (1978), the solar radiation force JF ra d on a surface 
area dA due to reflection and absorption at the solar panel 
is given by 


d F ra d = — P 


(1 — O s )s + 2(G S cos £ T- — Cd)n 


cos £dA, 
(A7) 


where C s and Cd are defined in Section Al, and n is the 


unit vector normal to the solar panel [fi = (0, 0, — 1)]. The 
solar radiation pressure P is 4.5 x 10' 6 Nm' 2 . We introduce 
the coefficients G and H\ 


<?=(!-GO, 


H = (1 + C s ) cos £ + -Ga, 


with which equation (A7) can be written as 


d-Gad = ~P 


G(— cos ip sin £ + 02) 
G(sin ip sin £ — ai) 
-H 


cos £dA. 


(AS) 


The solar radiation torque acting on the satellite is obtained 
by integrating the outer product of the force and the sepa¬ 
ration of the surface element from the CoG over the surface 


area of the solar panel. If (p, (p ) are polar coordinates defined 
relative to the centre of the solar panel, then the position 
vector r of the surface element d A is given by 


p cos cp + x 0 " 
p sin (p + yo , 
z 0 


(A9) 


where ( xo,yo,zo ) is the position of the centre of the solar 
panel relative to the CoG, measured in the SRS. The torque 
acting on the satellite is now given by 


N = 


V X d F rad, 


(A10) 


and the surface element by 


dA = pdpdcp. 


(All) 


In the integration, the terms containing tp in equation (A9) 
cancel, leaving only the effects of the offset of the CoG: 

0 z 0 1 r 

n COS Ip y- / a 1 

z 0 0 . . P, (A12) 

sin ip 

L —yo -Xo. L J 


N k,N 0 + — G sin 2£ 


where non-modulated part of the torque, Wo, is given by 


No = cos £ 


yoH — zoGai 
—xo H — zoGa 2 
G(x 0 ai +1/00:2) 


P, 


(A13) 


and the total force T = 7.16x 10 5 N. In the ISR the torques 
are given by 


No = cos£ 


yoH — z 0 Gai — xoHa 2 
—xo H — z 0 Ga 2 + yoHct3 
(G - H)(xoai + 1/002) 


P, 


(A14) 


while the product of sin 2^ and 0 1,2,3 are small enough to 
leave the second part in equation ( |A12| ) unchanged. 

The inertia tensor given in Section 0 defines an order 
of magnitude for the off-diagonal elements, based on the Ma- 
tra Marconi document, which provides estimates of 0.03 rad 
for oi and 02. The typical resulting solar radiation torques 
are at a level of 2 x 10~ 6 Nm in x and y, and 10~ 7 in z, 
equivalent to accelerations at a level of 1 to 0.1 milliarc- 
sec s' 2 . Scanning at a tilt of 10° adds a modulated torque 
with an amplitude of about 7.5 x 10 -6 Nm in x and y, and 
3.8 x 10'' Nm in z. In all cases, the effect of the solar radi¬ 
ation is small with respect to the beam sizes of the Planck 
detectors, and the only relevant contribution is the acceler¬ 
ation or deceleration of the spin rate, which will accumulate 
to give barely significant phase shifts over a one hour scan 
period. 


A5 Analytic model for the inertial rates 


In its slightly si mpli fied form, given by the first two terms 
only, equations (A4) can be solved to provide a model for 
the evolution of the inertial rates around the satellite axes. 
Assuming oo z t = ip, we find 

K, 0 


UJ X — 


Iyyf20J; 


— b x cosu> z f 


+wi cos(w 2 7t) — UJ 2 sin(u; z 7t), 

K,o 

IxxflOJi 


+ by sin ui z t 
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+ 


h L 


u >2 cos(to z ~/t) + uj\ sin(o> z 7 t) 
G sin 2£ 


, <o, 

W3 4- Tt—t — 


— (ya sin u> z t — xo cos a > z t)F, 
(A15) 

where Nf 0 , N' y 0 and N l n are the components of N o as de¬ 
fined in equation (A13), and lo z is the nominal scan velocity, 
and ( 07 , 07 , 07 ) are integration constants. Also 

_ 1 z 0 Gsin2£(2T yy /T zz - 1 ) 

* “ 2 


b v — — 


(in T lyy Izz)^z 

lz 0 G sin 2 £( 2 /' x /I' zz -\) 
2 


T, 


(Ixx T lyy J ZZ )o7 


T, 


(A16) 


where F, as defined in equation (A14), is the total solar ra¬ 
diation force, and u > z 7 is the nutation frequency, equivalent 
to a nutation period of just under 3 minutes for 7 = 0.35 
and 07 equal to 1 rpm. The satellite will undergo nutation 
damping after each repositioning of the spin axis. This re¬ 
duces the rotation rates around the x and y axes to a level 
low enough not to cause significant disturbances on the mea¬ 
surements. If we assumed that, as a result of the nutation 
damping, the rotation rates at the start of a scan are equal 
to ( 07 , 0 , 07 ,, 0 , 07 , 0 ), then ( 07 , 07 , 07 ) are given by 

K, 0 


07 


UJ2 


Wx,0 + 


h 

7 L 


I, 

(-Oy ,0 


yyU} z f2 

Nf 0 


+ bx 


Ixx^zfl 


G sin 2f 

07 = 07 , 0 — Xo T, -F 1 

2a izizz 


(AIT) 


where the reference time t = 0 is the first instance of ip = 0 
after the end of the nutation damping process. 


A6 Analytic model for the error angles 

The differences between the nominal and actual orientations 
of the IRS axes are described by a set of angles generally re¬ 
ferred to as the Tait-Bryan angles (T z ,T y ,T x ). The relations 
between increments in the Tait-Bryan angles and the iner¬ 
tial rates are described by the equations give below, where it 
is assumed that the inertial rates resulting from the nominal 
scanning law are given by 07 , and the actual inertial rates 
by (07 , a J y , uj z ). 

The inertial rate for the nominal scanning law needs to 
be transformed to inertial rates around the actual satellite 
axes: 


(AJx,i 


" sin Tx sin T z — cos T x cos T z sin T y " 


— Uz 

sin Tx cos T z + cos T x sin T z sin T y 



_ cos Tx cos Ty 


(A 18) 


The differences between the actual rates and the projected 
nominal rates are transformed to corrections to the Tait- 
Bryan angles: 

cos T z / cos T v — sin T z / cos T y 
sin T z cos T z 


T = 


— cos T z tan T v 


sin T z tan T v 


Slj, 


where 


= 


(A19) 


(A20) 
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The first set of equations can be approximated by 


(A)x,i 


'~Ty' 


w 0 ) z 

T x 

_W z ,i_ 
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(A21) 


Similarly, the second set of equations becomes 

T x ss (wi + T y <jj z ), 

Ty (O Jy TeO^z), 


T z 


(a 7 07 ), 


(A22) 


where the relations for (uj x ,u y ) given in equation (4.15) 
can be use d to obtain an approximate solution for (T x ,T y ). 
Equation (A15) is written as 


Mx 

UJy 


= —a x — b x cosa ) z t + 07 cos(o77f) — 07 sin(o77i), 
= a y +by sin u> z t 


+ 


[a >2 cos(a77 t) -I- 07 sin(a77t)]. 


(A23) 


Similarly, the angles ( T x ,Ty ) are represented by 
T x = A x + B x cos ui z t + C x sin a i z t + tF cos a j z t 
+D X cos Q z jt + E x sin 077 !, 

= A y + By cos ui z t — Cy sin Ld z t — tF sin a J z t 
+D y cos uiz-yt + E y sin 077 !. 


T, 


(A24) 


Substituting equations (A23) and (A24) in equation (A22), 
we find 


D x = 


E x — 


F = 


A X Ojy !LO Z 

C x = Ci 
B x =C 2 
^2 ( 7 //l - 7 ) 

07(1 - 7 2 ) 
^ 1 ( 7//1 ~ 7) 

u>*(l - 7 2 ) 

— (5x + by) 


A y — n x j 07 , 

D _ (bx — by) , ^ 

t>y - I 01 , 

2iUz 

Cy = C 2 , 

_ 07 (/ 2 - 1) 

" 07(1 — 7 2 ) ’ 

_ Q7 (/ 2 - 1) 

V 0,(1- 7 2 )’ 


(A25) 


where C\ is a constant, depending on the choice of the ref¬ 
erence position of the satellite’s 2 -axis. If ( T x fi,Ty,o ) are the 
Tait-Bryan angles at time zero, then for 5 = 0 the following 
relations are obtained: 


Tx. 0 = 


T y ,o = 


ay , ^ , 07 ( 7 /fi — 7 ) 
07 w z (l — 7 2 ) 

Ox n 07 (/ 2 - 1) 

07 07(1 — 7 2 ) 


(A26) 


where a x and a y are relatively small with respect to the last 
terms. It is possible to choose the reference position for the 
Tait-Bryan angles such that Ci and C 2 are equal to zero. 
This reduces the description of the motion of the satellite 
2 -axis to a simple ellipse in the (T x ,T y ) plane. 

The coefficient for F is a potential source of instabil¬ 
ity. An ev aluation of the expected value of F, using equa¬ 
tion (A16), gives 


_ zqG sin 2£ r 

2 0 , 


A1 


(A27) 


shows 


Using the input parameters specified in Section 
that the effect is quite small. From the values given we esti¬ 
mate a maximum drift of 166 x sin 2^ arcsec over a 1 hour 
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Figure Al. A typical evolution pattern for inertial rates and error angles for £ = 10 degrees. The two graphs in the upper right and 
lower left corners show the almost perfect relations between rates and error angles. 


period. With £ having a maximum value of 10 degrees, this 
is noticeable at the resolutions used by the Planck instru¬ 
ments. 

It can be verified that for £ = 0 and B = C\ = C 2 = 0, 
and a fixed value of uj z , linear relations exist between T x and 
uj y and between T y and uj x . The ratios are 


R _ T * = 

UJy 

n _ Ty — CLx/^z 


1 1-/1 
Wz 1 - 7 2 ’ 

-1 1-/2 
uJz 1 - 7 2 ‘ 


(A28) 


These relations, shown in Fig. Al, get only slightly disturbed 
by non-zero values for the solar aspect angle. They show 
that the satellite dynamics is dominated by the effects of 
the residual velocities around the x and y axes. 


A7 Numerical integration of the attitude model 

Numerical integration of the Euler equation and the result¬ 
ing rates can be performed without any of the approxima¬ 
tions made in the analytic model; this provides a useful cross 
check on the analytic results. The integration over the Eu¬ 
ler equation requires a set of inertial rates at the starting 
time. These have been chosen such that they represent the 
nominal rates: zero in uj x and u y , 6 degrees per second in 
ui z . At time-steps of 0.01 s the torque acting on the satel¬ 
lite is developed, as well as the crossproducts of the rates 
l _1 cu x I u). This integration provides a record of the inertial 
rate development. 


Similarly, the integration over the Tait-Bryan angles re¬ 
quires an ass umpt ion about the starting values. For T x and 
T y equation (A2f) is used to calculate the starting values, 
with Ci = Ci = 0, while the T, starting value is set to 




A 7.1 Comparison with the analytic results 


The results of the analytic approximation are compared with 
the numerical integrations in Fig. AS. These comparisons 
show that over a one-hour period the analytic approxima¬ 
tion adds noise at a level of a few arcsec to the more accu¬ 
rate numerical integration. The main component of the noise 
originates from ignoring the I' xy term in the development of 
equation (A4). Including this term affects the amplitudes 


wi and 0 J 2 , and adds to these a small time-dependent term, 
which can, if necessary, also be estimated analytically. How¬ 
ever, these corrections appear to be insignificant for Planck. 
The positional error in the analytic model is of the order of 
a few arcsec, which again is insignificant for Planck. 

Given this level of agreement between the analytic and 
numerical models, it will be sufficient for the Planck simu¬ 
lations to use the analytic model, and to simulate only the 
initial conditions for a scan. In the same way, the condi¬ 
tions during re-positioning of the 2 -axis can be simulated. 
The simulations, as well as the analytic model, show that a 
solar aspect angle of 10° or less has very little influence on 
the satellite’s attitude, unless this is examined at accuracies 
much higher than those required for Planck. 
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time (s) 

Figure A2. The differences between rates and error angles around the x and y axes as obtained from analytic modelling and from 
numerical integration. Top two graphs: a one-hour scan with £ = 0; bottom two graphs: £ = 10 degrees. The two curves in each graph 
show the variations for the x and y axes. 


APPENDIX B: PROJECTIONS ON THE 
REFERENCE CIRCLE 

In this appendix we describe the transformations for posi¬ 
tions and rates from a system in which the spin axis de¬ 
scribes a nutation ellipse in the satellite reference frame to 
a reference circle with its axis at the centre of this ellipse. 

For £ = 0°, the nutation ellipse is described in the satel¬ 
lite reference system by 

T x « A x + D x cos ujz'yt + E x sin uj^t 

Ty Ri Ay + Dy COS U) Z ^t + Ey sin ui z -yt , (Bl) 

© 2001 RAS, MNRAS 000, 1 |(] 


with the c oefficients as defined in equation (A2i). In equa¬ 
tion (A17) it can in addition be assumed that (co Xl o, ui y ,o) 
are by far the dominant coefficients in (u}i,u> 2 )- This also 
implies that the coefficients (A x ,A y ) will be relatively small 
with respect to ( 101 , 012 ), so that we are left with 


T x 


£ 

L 0 x ncos(u) z 'ft) — —LOynSm(Co z -yt), 

7 

'y 

Uyfi cos(u)z7f) + — L 0 X}0 sin(u> z -ft), 

J 1 

Ri [ojyfi cos(o)z7f) + j-oj Xi o sin(cu z 7t)], 
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£ 

Ty « R 2 \uj x , 0 cos(u> z -yt) - —w. V: o sin(±z 7 t)], (B2) 

where 7?i and R 2 are functions of /i, f 2 and ui z , as defined 


in equation (A28), and 7 = \Jfif 2 as defined in Section A3 


By choosing a suitable reference time to equation (B2) can 
be written as 


(jJ x r- 

a uj x ,o cos [w z 7 (t — to)], 


U)y ? 

a j-w®, osin[w* 7 (£-£ 0 )], 


T x ? 

a /Zi j-w*, 0 sin[w*7(t —to)], 


Ty P 

a R 2 U> x ,ocos[ui z y(t — to)], 

(B3) 


The mean spin axis is assumed to be a fixed point in space, 
in the centre of the ellipse described by ( T x ,T y ), and the 
spin velocity is assumed to be the nominal spin velocity, 
t j z = ui z . The offset angles can be defined with respect to a 
fixed reference system relative to the mean spin axis: 

T' x = T x cos ip — T y sin tp, 

Ty = T x sin tp + T y cos tp. (B4) 


The angle tp is measured on the ring from its crossing of the 
reference circle (see Fig. []]). A point on the ring is given by 
the direction cosines 


Oi = 


cos tp sin a 
sin tp sin a 
cos a 


(B5) 


where a is the opening angle for the reference circle. The 
offset angles (T' x , Ty) create small offsets in the actual angles 
(tp, a) for a measurement. We ignore the very small rotation 
over T z . The rotations over T x and Ty distort equation (|B5|) 


Oi = 


cos tp sin a + Ty cos a 
sin tp sin a — T' x cos a 
cos a + T X sin a sin tp — Ty sin a sin tp 


(B 6 ) 


which can also be expressed as 


Oi = 


cos tp sin a + Aa cos a cos tp — Atp sin a sin tp 
sin tp sin a + Aa cos a sin tp ± Atp sin a cos tp 
cos a — Aa sin a 


. (B7) 


From these relations the expressions for (Atp, Aa) in terms 
of (T' x ,Ty) follow: 

Atp — (—T x cosip — Ty simp)/ tana, 

Aa = —T x sin tp + T' y cos tp. (B 8 ) 


Substituting equation (B4) we find 


Atp = —T x / tana, 

Aa = +T y . (B9) 


With an opening angle a of 80 to 85 degrees, the effect of 
the spin-axis wobble on the scan phases is very small. The 
maximum amplitude is expected to be less than 0.12 arcmin. 
This may be visible, though, in the star mapper data. The 
effect of T x on the orientation of the focal plane is even 
smaller. 

The a variations can be relevant for the high-frequency 
detectors. A change in position of the beam by ±0.75 arcmin 
for a point source close to the steepest slope of the beam 
will create significant variations in the signal over a TOP. If 
the wobble parameters can be reconstructed from the star 


mapper data, then some of this variation may be accounted 
for in the reductions. 

This paper has been produced using the Royal Astronomical 
Society/Blackwell Science IAT^K style file. 


© 2001 RAS, MNRAS 000, id 












